




#ifndef _COMPLEX_H_
#define _COMPLEX_H_




#include "stl_type.h"



// 这里实现复数的运算，FFT，解一元n次方程（可以求解复数根）


/// <summary>
/// C语言接口的向量结构体
/// </summary>
typedef struct _complex__
{
	real_t real;
	real_t imag; //imaginary虚部
}complex;

/// <summary>
/// 打印复数
/// </summary>
/// <param name="c"></param>
void complex_print(complex c);

/// <summary>
/// 判断复数是否相等。复数没法直接比大小。
/// </summary>
/// <param name="A">复数</param>
/// <param name="B">复数</param>
/// <param name="precision">精度</param>
/// <returns>A与B的差在精度范围内，返回true，否则返回false</returns>
bool complex_isEqual(complex A, complex B, real_t precision);

/// <summary>
/// 判断是否是零。 
/// </summary>
/// <param name="A">复数</param>
/// <param name="precision">精度</param>
/// <returns>A与0的差在精度范围内，返回true，否则返回false</returns>
bool complex_isZero(complex A, real_t precision);

/// <summary>
/// 计算复数模长
/// </summary>
/// <param name="A">复数</param>
/// <returns>复数的模长</returns>
real_t complex_mod(complex A);
/// <summary>
/// 
/// </summary>
/// <param name="A">复数</param>
/// <returns>返回复数的模的平方</returns>
real_t complex_modSquare(complex A);
/// <summary>
/// 
/// </summary>
/// <param name="A">复数</param>
/// <returns>共轭复数</returns>
complex complex_conj(complex A);
/// <summary>
/// 复数加法
/// </summary>
/// <param name="a">复数</param>
/// <param name="b">复数</param>
/// <returns>两个复数的和 a+b</returns>
complex complex_add(complex a, complex b);
/// <summary>
/// 复数减法
/// </summary>
/// <param name="a">复数</param>
/// <param name="b">复数</param>
/// <returns>两个复数的差 a-b</returns>
complex complex_sub(complex a, complex b);
/// <summary>
/// 复数乘法
/// </summary>
/// <param name="a">复数</param>
/// <param name="b">复数</param>
/// <returns>两个复数的积 a*b</returns>
complex complex_mul(complex a, complex b);

/// <summary>
/// 复数乘实数
/// </summary>
/// <param name="a">复数</param>
/// <param name="b">实数</param>
/// <returns> a * b </returns>
complex complex_mul_real(complex a, real_t b);


/// <summary>
/// 复数除法
/// </summary>
/// <param name="a">复数</param>
/// <param name="b">复数</param>
/// <returns>两个复数的商 a/b</returns>
complex complex_div(complex a, complex b);

/// <summary>
/// 复数除实数
/// </summary>
/// <param name="a">复数</param>
/// <param name="b">实数</param>
/// <returns>复数的商 a/b</returns>
complex complex_div_real(complex a, real_t b);

//根据time_domain计算frequency_domain,w 是旋转因子缓存，长度是length的一半。

/// <summary>
/// 根据time_domain计算frequency_domain,w 是旋转因子缓存，长度是length的一半。
/// </summary>
/// <param name="time_domain">时域</param>
/// <param name="frequency_domain">频域</param>
/// <param name="w">旋转因子，中间计算量的缓存，长度是length的一半</param>
/// <param name="log2_length">长度的2为底的对数。</param>
void fft(complex* time_domain, complex* frequency_domain, complex* w, int log2_length);

/// <summary>
/// 根据frequency_domain计算time_domain
/// </summary>
/// <param name="time_domain">时域</param>
/// <param name="frequency_domain">频域</param>
/// <param name="w">旋转因子，中间计算量的缓存，长度是length的一半</param>
/// <param name="log2_length">长度的2为底的对数。</param>
void ifft(complex* time_domain, complex* frequency_domain, complex* w, int log2_length);


/// <summary>
/// 测试傅里叶变换
/// </summary>
void test_FT(void);


typedef struct _wiener__
{
	/// <summary>
	///  长度的2为底的对数。
	/// </summary>
	int log2_length;
	/// <summary>
	/// 信号数组的长度
	/// </summary>
	int length;
	/// <summary>
	/// 信号样本
	/// </summary>
	complex* signal;
	/// <summary>
	/// 信号频谱
	/// </summary>
	complex* signal_frequency;
	/// <summary>
	/// 噪声样本
	/// </summary>
	complex* noise;
	/// <summary>
	/// 噪声频谱
	/// </summary>
	complex* noise_frequency;

	/// <summary>
	/// 观测样本
	/// </summary>
	complex* observe;
	/// <summary>
	/// 观测频谱
	/// </summary>
	complex* observe_frequency;
	/// <summary>
	/// 噪声模长
	/// </summary>
	real_t* noise_mod;

	/// <summary>
	/// 傅里叶变换的缓存数组。长度是length / 2
	/// </summary>
	complex* ftbuf;

}wiener;


/// <summary>
/// 分配维纳滤波器内存
/// </summary>
/// <param name="w">维纳滤波器</param>
/// <param name="ln2len">滤波器处理的信号长度</param>
void wiener_allocate(wiener* w, int ln2len);
/// <summary>
/// 释放维纳滤波器内存
/// </summary>
/// <param name="w">维纳滤波器</param>
void wiener_release(wiener* w);

/// <summary>
/// 用噪声样本初始化维纳滤波器
/// </summary>
/// <param name="w">维纳滤波器</param>
/// <param name="noise">噪声数组</param>
/// <param name="len">数组长度</param>
void wiener_init(wiener w, real_t* noise, int len);


/// <summary>
/// 对维纳滤波器进行滤波。
/// 先在外部更新observe数组。
/// 然后滤波器先计算observe的频谱，减去噪声频谱后得到signal频谱，再根据signal频谱计算signal。
/// 滤波后的数据会放到signal数组。
/// </summary>
/// <param name="w">维纳滤波器</param>
void wiener_filter(wiener w);

//测试维纳滤波。 只要将原始数据，噪声频谱，滤波数据都打印出来，excel画图即可。
void test_weina(void);


//这里用 两个主频率的噪声，加上sigmoid函数。进行测试。
void test_weina2(void);








#endif







//

